package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.TP3;
import static scu.maqiang.numeric.Constants.TW3;
import static scu.maqiang.numeric.MVO.*;

import java.util.Arrays;

import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.Constants;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.ParamCheck;
import scu.maqiang.numeric.SRMatrix;

public class FES2T32 extends FES {
	public FES2T32(Mesh2T3 mesh2D) {
		super(mesh2D, 2);
		eleStiff = FES2T32::elementElasticityStiff;
		eleMass = FES2T32::elementElasticMass;
		eleSource = FES2T32::elementVolumeForce;
		eleFlux = FES2T32::elemenSurfaceForce;
	}

	/**
	 *
	 * @param matCoef
	 * @param tp
	 * @param A
	 */
	public void assembleElasticStiff(double[][] matCoef, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(matCoef, FES2T32::elementElasticityStiff, tp, A);
	}

	/**
	 *
	 * @param coefFunc
	 * @param param
	 * @param tp
	 * @param A
	 */
	public void assembleElasticStiff(MatrixFunc coefFunc, double[][][] param, BVPType tp, SRMatrix A) {
		assembleGlobalMatrix(coefFunc, param, FES2T32::elementElasticityStiff, tp, A);
	}

	/**
	 *
	 * @param xy
	 * @param ENu
	 * @param tp
	 * @param K
	 */
	public static void elementElasticityStiff(int[] ele, double[][] xy, double[][] ENu, BVPType tp, double[][] K) {
		double[] x = new double[] { xy[0][0], xy[1][0], xy[2][0] };
		double[] E = ENu[0];
		double[] Nu = ENu[1];
		double[][] grad = new double[3][2];
		double area = CG.gradient2DT3(xy, grad);
		//System.out.println(MVO.toString(xy));
		//System.out.println(MVO.toString(grad));
		double mE = average(E);
		double mNu = average(Nu);
		double[][] B = null;
		double[][] D = null;
		double r0 = average(x);
		MVO.fill(K, 0.0);
		switch (tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
			D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
			B = new double[3][6];
//			B[0][0] = B[2][1] = grad[0][0];
//			B[0][2] = B[2][3] = grad[1][0];
//			B[0][4] = B[2][5] = grad[2][0];
//			B[1][1] = B[2][0] = grad[0][1];
//			B[1][3] = B[2][2] = grad[1][1];
//			B[1][5] = B[2][4] = grad[2][1];
			for(int ib = 0; ib < 3; ib++) {
				B[0][2 * ib] = B[2][2 * ib + 1] = grad[ib][0];
				B[1][2 * ib + 1] = B[2][2 * ib] = grad[ib][1];
			}
			mul(matmul(transpose(B), matmul(D, B)), area, K);
			break;
		case CONSTITUTE_MATRIX_COMMON:
			B = new double[3][6];
//			B[0][0] = B[2][1] = grad[0][0];
//			B[0][2] = B[2][3] = grad[1][0];
//			B[0][4] = B[2][5] = grad[2][0];
//			B[1][1] = B[2][0] = grad[0][1];
//			B[1][3] = B[2][2] = grad[1][1];
//			B[1][5] = B[2][4] = grad[2][1];
			for(int ib = 0; ib < 3; ib++) {
				B[0][2 * ib] = B[2][2 * ib + 1] = grad[ib][0];
				B[1][2 * ib + 1] = B[2][2 * ib] = grad[ib][1];
			}
			mul(matmul(transpose(B), matmul(ENu, B)), area, K);
			break;
		case AXISYMMETRIC:
			D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
			B = new double[4][6];
			area = area * r0;
//			B[0][0] = B[3][1] = grad[0][0];
//			B[0][2] = B[3][3] = grad[1][0];
//			B[0][4] = B[3][5] = grad[2][0];
//			B[1][0] = B[1][2] = B[1][4] = 1.0 / (3 * r0);
//			B[2][1] = B[3][0] = grad[0][1];
//			B[2][3] = B[3][2] = grad[1][1];
//			B[2][5] = B[3][4] = grad[2][1];
			for(int ib = 0; ib < 3; ib++) {
				B[0][2 * ib] = B[3][2 * ib + 1] = grad[ib][0];
				B[1][2 * ib + 1] = B[3][2 * ib] = grad[ib][1];
				B[2][2 * ib] = 1.0 / (3.0 * r0);
			}
			mul(matmul(transpose(B), matmul(D, B)), area, K);
			break;
		case CONSTITUTE_MATRIX_AXIS:
			B = new double[4][6];
			area = area * r0;
//			B[0][0] = B[3][1] = grad[0][0];
//			B[0][2] = B[3][3] = grad[1][0];
//			B[0][4] = B[3][5] = grad[2][0];
//			B[1][0] = B[1][2] = B[1][4] = 1.0 / (3 * r0);
//			B[2][1] = B[3][0] = grad[0][1];
//			B[2][3] = B[3][2] = grad[1][1];
//			B[2][5] = B[3][4] = grad[2][1];
			for(int ib = 0; ib < 3; ib++) {
				B[0][2 * ib] = B[3][2 * ib + 1] = grad[ib][0];
				B[1][2 * ib + 1] = B[3][2 * ib] = grad[ib][1];
				B[2][2 * ib] = 1.0 / (3.0 * r0);
			}
			mul(matmul(transpose(B), matmul(ENu, B)), area, K);
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public static void elementElasticMass(int[] ele, double[][] xy, double[][] coef, BVPType tp, double[][] K) {
		double area = CG.area2DT3(xy);
		double[] ccoef = coef[0].clone();
		MVO.fill(K, 0.0);
		double mcoef = 0.0;
		if (tp == BVPType.CONSIST_MASS) {
			for(int i = 0; i < 3; i++) {
				ccoef[i] *= area / 60.0;
			}
			mcoef = MVO.sum(ccoef);
			K[0][0] = K[1][1] = 2 * (2 * ccoef[0] + mcoef);
			K[2][2] = K[3][3] = 2 * (2 * ccoef[1] + mcoef);
			K[4][4] = K[5][5] = 2 * (2 * ccoef[2] + mcoef);
			K[0][2] = K[2][0] = K[1][3] = K[3][1] = (ccoef[0] + ccoef[1] + mcoef);
			K[0][4] = K[4][0] = K[1][5] = K[5][1] = (ccoef[0] + ccoef[2] + mcoef);
			K[2][4] = K[4][2] = K[3][5] = K[5][3] = (ccoef[1] + ccoef[2] + mcoef);
			
			//System.out.println(MVO.toString(elemMatrix));
		} else if (tp == BVPType.LUMP_MASS) {
			for(int i = 0; i < 3; i++) {
				ccoef[i] *= area / 12.0;
			}
			mcoef = MVO.sum(ccoef);
			K[0][0] = K[1][1] = ccoef[0] + mcoef;
			K[2][2] = K[3][3] = ccoef[1] + mcoef;
			K[4][4] = K[5][5] = ccoef[2] + mcoef;
		} else if (tp == BVPType.AXIS_CONSIST_MASS) {
			for(int i = 0; i < 3; i++) {
				ccoef[i] *= xy[i][0] * area / 60.0;
			}
			mcoef = MVO.sum(ccoef);
			K[0][0] = K[1][1] = 2 * (2 * ccoef[0] + mcoef);
			K[2][2] = K[3][3] = 2 * (2 * ccoef[1] + mcoef);
			K[4][4] = K[5][5] = 2 * (2 * ccoef[2] + mcoef);
			K[0][2] = K[2][0] = K[1][3] = K[3][1] = (ccoef[0] + ccoef[1] + mcoef);
			K[0][4] = K[4][0] = K[1][5] = K[5][1] = (ccoef[0] + ccoef[2] + mcoef);
			K[2][4] = K[4][2] = K[3][5] = K[5][3] = (ccoef[1] + ccoef[2] + mcoef);
		} else if (tp == BVPType.AXIS_LUMP_MASS) {
			for(int i = 0; i < 3; i++) {
				ccoef[i] *= xy[i][0] * area / 12.0;
			}
			mcoef = MVO.sum(ccoef);
			K[0][0] = K[1][1] = ccoef[0] + mcoef;
			K[2][2] = K[3][3] = ccoef[1] + mcoef;
			K[4][4] = K[5][5] = ccoef[2] + mcoef;
		} else {
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public void assembleVolumeForce(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS) {
		assembleGlobalVector(coefFunc, param, FES2T32::elementVolumeForce, tp, RHS);
	}

	public void assembleVolumeForce(double[] constCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(constCoef, FES2T32::elementVolumeForce, tp, RHS);
	}

	public static void elementVolumeForce(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double area = CG.area2DT3(xy);
		double[] fu = f[0];
		double[] fv = f[1];
		double[] x = new double[] { xy[0][0], xy[1][0], xy[2][0] };
		double sumFu = MVO.sum(f[0]);
		double sumFv = MVO.sum(f[1]);
		Arrays.fill(eleF, 0.0);
		switch (tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
		case COMMON:
		case CONSTITUTE_MATRIX_COMMON:
			for (int i = 0; i < 3; i++) {
				eleF[2 * i] = (fu[i] + sumFu) * area / 12.0;
				eleF[2 * i + 1] = (fv[i] + sumFv) * area / 12.0;
			}
			break;
		case CONSTITUTE_MATRIX_AXIS:
		case AXISYMMETRIC:
			for (int j = 0; j < 4; j++) {
				for (int i = 0; i < 3; i++) {
					eleF[2 * i] += area * TW3[j] * MVO.dot_product(TP3[j], fu) * MVO.dot_product(TP3[j], x) * TP3[j][i];
					eleF[2 * i + 1] += area * TW3[j] * MVO.dot_product(TP3[j], fv) * MVO.dot_product(TP3[j], x)
							* TP3[j][i];
				}
			}
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public void assembleSurfaceForce(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS,
			int[] bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES2T32::elemenSurfaceForce, tp, RHS, bdLabel);
	}

	public void assembleSurfaceForce(double[] constCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(constCoef, FES2T32::elemenSurfaceForce, tp, RHS, bdLabel);
	}
	
	public void assembleSurfaceForce(double[][] vecCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(vecCoef, FES2T32::elemenSurfaceForce, tp, RHS, bdLabel);
	}

	public static void elemenSurfaceForce(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double length = CG.length2DL2(xy);
		double[] gu = f[0];
		double[] gv = f[1];
		double sumGU = MVO.sum(gu);
		double sumGV = MVO.sum(gv);
		switch (tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
		case COMMON:
		case CONSTITUTE_MATRIX_COMMON:
			length /= 6.0;
			for (int i = 0; i < 2; i++) {
				eleF[2 * i] = (sumGU + gu[i]) * length;
				eleF[2 * i + 1] = (sumGV + gv[i]) * length;
			}
			break;
		case AXISYMMETRIC:
		case CONSTITUTE_MATRIX_AXIS:
			length /= 12.0;
			double[] x = MVO.col(xy, 0);
			double sumX = MVO.sum(x);
			for (int i = 0; i < 2; i++) {
				eleF[2 * i] = (sumGU * sumX + 2 * gu[i] * x[i]) * length;
				eleF[2 * i + 1] = (sumGV * sumX + 2 * gv[i] * x[i]) * length;
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
		
	}

	public void assemblePressure(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS,
			int[] bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES2T32::elementPressure, tp, RHS, bdLabel);
	}

	public void assemblePressure(double[] constCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(constCoef, FES2T32::elementPressure, tp, RHS, bdLabel);
	}
	
	public void assemblePressure(double[][] vecCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(vecCoef, FES2T32::elementPressure, tp, RHS, bdLabel);
	}
	
	public static void elementPressure(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double[] pressure = f[0];
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[][] ff = new double[2][2];
		double length = Math.sqrt((x[0] - x[1])*(x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]));
		double nx = (y[1] - y[0])/length;
		double ny = (x[0] - x[1])/length;
		ff[0][0] = -pressure[0] * nx;
		ff[0][1] = -pressure[1] * nx;
		ff[1][0] = -pressure[0] * ny;
		ff[1][1] = -pressure[1] * ny;
		elemenSurfaceForce(ele, xy, ff, tp, eleF);
	}	

	public void assembleHeatLoad(ScalarFunc[] coef, double[] param, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, param, vecCoef, FES2T32::elementHeatLoad, tp, RHS);
	}

	public void assembleHeatLoad(double[] coef, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, vecCoef, FES2T32::elementHeatLoad, tp, RHS);
	}

	public static void elementHeatLoad(int[] ele, double[][] coord, double[][] coef, BVPType tp, double[] eleLoad) {
		double[] beta = null;
		double[] temperature = null;
		double[] x = MVO.col(coord, 0);
		double[][] grad = new double[3][2];
		double[] dxN = new double[3];
		double[] dyN = new double[3];
		double[] N = new double[3];
		double area = CG.gradient2DT3(coord, grad);
		Arrays.fill(eleLoad, 0.0);
		switch (tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
		case COMMON:
			beta = coef[0];
			temperature = coef[1];
			double c = (MVO.sum(beta) * MVO.sum(temperature) + MVO.dot_product(beta, temperature)) * area / 12.0;
			for (int i = 0; i < 3; i++) {
				eleLoad[2 * i] = c * grad[i][0];
				eleLoad[2 * i + 1] = c * grad[i][1];
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			double[] beta1 = coef[0];
			double[] beta2 = coef[1];
			temperature = coef[2];
			double cx = (MVO.sum(beta1) * MVO.sum(temperature) + MVO.dot_product(beta1, temperature)) * area / 12.0;
			double cy = (MVO.sum(beta2) * MVO.sum(temperature) + MVO.dot_product(beta2, temperature)) * area / 12.0;
			for (int i = 0; i < 3; i++) {
				eleLoad[2 * i] = cx * grad[i][0];
				eleLoad[2 * i + 1] = cy * grad[i][1];
			}
			break;
		case AXISYMMETRIC:
			beta = coef[0];
			temperature = coef[1];
			for (int i = 0; i < 4; i++) {
				double[] L = Constants.TP3[i];
				double wxarea = MVO.dot_product(beta, L) * MVO.dot_product(temperature, L) * area * TW3[i];
				for (int j = 0; j < 3; j++) {
					eleLoad[2 * j] += (grad[j][0] * MVO.dot_product(x, L) + L[j]) * wxarea;
					eleLoad[2 * j + 1] += grad[j][1] * MVO.dot_product(x, L) * wxarea;
				}
			}
			break;
		case CONSTITUTE_MATRIX_AXIS:
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
//		System.out.println(Arrays.toString(eleLoad));
//		System.exit(0);
	}

	public void computeNodalStrain(double[] sol, double[][] strain) {
		int len = sol.length;
		ParamCheck.checkEqual(len, this.ndof);
		int n = len / 2;
		ParamCheck.checkEqual(n, strain[0].length);
		ParamCheck.checkEqual(n, strain[1].length);
		ParamCheck.checkEqual(n, strain[2].length);
		double[] angle = new double[3];
		double[][] coord = new double[3][2];
		double[][] grad = new double[3][2];
		double[] angleCount = new double[this.GetNdof()];
		double dxU, dyU, dxV, dyV;
		double[] eleU = new double[3];
		double[] eleV = new double[3];
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			CG.gradient2DT3(coord, grad);
			CG.angle2DT3(coord, angle);
			for (int j = 0; j < 3; j++) {
				eleU[j] = sol[2 * ele[j]];
				eleV[j] = sol[2 * ele[j] + 1];
			}
			double[] dx = MVO.col(grad, 0);
			double[] dy = MVO.col(grad, 1);
			dxU = MVO.dot_product(eleU, dx);
			dyU = MVO.dot_product(eleU, dy);
			dxV = MVO.dot_product(eleV, dx);
			dyV = MVO.dot_product(eleV, dy);
			for (int j = 0; j < 3; j++) {
				strain[0][ele[j]] += dxU * angle[j];
				strain[1][ele[j]] += dyV * angle[j];
				strain[2][ele[j]] += (dxV + dyU) * angle[j];
				angleCount[ele[j]] += angle[j];
			}
		}
		for (int i = 0; i < n; i++) {
			strain[0][i] /= angleCount[i];
			strain[1][i] /= angleCount[i];
			strain[2][i] /= angleCount[i];
		}
	}

	public void computeElementStrain(double[] sol, double[][] strain) {
		int len = sol.length;
		ParamCheck.checkEqual(len, this.ndof);
		int n = len / 2;
		int nt = mesh.getNt();
		ParamCheck.checkEqual(nt, strain[0].length);
		ParamCheck.checkEqual(nt, strain[1].length);
		ParamCheck.checkEqual(nt, strain[2].length);
		double[] angle = new double[3];
		double[][] coord = new double[3][2];
		double[][] grad = new double[3][2];
		double dxU, dyU, dxV, dyV;
		double[] eleU = new double[3];
		double[] eleV = new double[3];
		for (int i = 0; i < nt; i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			CG.gradient2DT3(coord, grad);
			CG.angle2DT3(coord, angle);
			for (int j = 0; j < 3; j++) {
				eleU[j] = sol[2 * ele[j]];
				eleV[j] = sol[2 * ele[j] + 1];
			}
			double[] dx = MVO.col(grad, 0);
			double[] dy = MVO.col(grad, 1);
			dxU = MVO.dot_product(eleU, dx);
			dyU = MVO.dot_product(eleU, dy);
			dxV = MVO.dot_product(eleV, dx);
			dyV = MVO.dot_product(eleV, dy);
			strain[0][i] = dxU;
			strain[1][i] = dyV;
			strain[2][i] = dxV + dyU;
		}
	}

	public void computeNodalStress(double[] sol, double[][] stress) {

	}

	public void computeElementStress(double[] sol, double[][] stress) {

	}

	public void computeNodalStrainStressVonMisesStress(BVPType tp, double[] coef, double[] sol, double[][] strain,
			double[][] stress, double[] misesStress) {
		int len = sol.length;
		ParamCheck.checkEqual(len, this.ndof);
		int n = len / 2;
		ParamCheck.checkEqual(n, strain[0].length);
		ParamCheck.checkEqual(n, strain[1].length);
		ParamCheck.checkEqual(n, strain[2].length);
		double[] angle = new double[3];
		double[][] coord = new double[3][2];
		double[][] grad = new double[3][2];
		double[] angleCount = new double[this.GetNdof()];
		double dxU, dyU, dxV, dyV;
		double[] eleU = new double[3];
		double[] eleV = new double[3];
		double E = coef[0];
		double Nu = coef[1];
		double mu = E / (2 * (1 + Nu));
		double lambda;
		double ex, ey, exy, ez, tx, ty, tz, txy, ttheta;
		double t1, t2, t3;
		double meanStress;
		if (tp == BVPType.PLANESTRESS) {
			lambda = E * Nu / (1 - Nu * Nu);
		} else {
			lambda = E * Nu / ((1 - 2 * Nu) * (1 + Nu));
		}
		for (int i = 0; i < mesh.getNt(); i++) {
			int[] ele = mesh.getElement(i);
			mesh.getCoordPerElement(ele, coord);
			CG.gradient2DT3(coord, grad);
			CG.angle2DT3(coord, angle);
			for (int j = 0; j < 3; j++) {
				eleU[j] = sol[2 * ele[j]];
				eleV[j] = sol[2 * ele[j] + 1];
			}
			double[] x = MVO.col(coord, 0);
			double[] dx = MVO.col(grad, 0);
			double[] dy = MVO.col(grad, 1);
			dxU = MVO.dot_product(eleU, dx);
			dyU = MVO.dot_product(eleU, dy);
			dxV = MVO.dot_product(eleV, dx);
			dyV = MVO.dot_product(eleV, dy);
			// for(int j = 0; j < 3; j++) {
			// strain[0][ele.idx[j]] += dxU * angle[j];
			// strain[1][ele.idx[j]] += dyV * angle[j];
			// strain[2][ele.idx[j]] += (dxV + dyU) * angle[j];
			// if (BVPType.AXISYMMETRIC == tp) {
			// strain[3][ele.idx[j]] += eleU[j] / x[j] * angle[j];
			// stress[0][ele.idx[j]] += (lambda + 2 * mu) * strain[0][ele.idx[j]] + lambda *
			// (strain[1][ele.idx[j]] + strain[3][ele.idx[j]]);
			// stress[1][ele.idx[j]] += lambda * (strain[0][ele.idx[j]] +
			// strain[3][ele.idx[j]]) + (lambda + 2 * mu) * strain[1][ele.idx[j]];
			// stress[3][ele.idx[j]] += (lambda + 2 * mu) * strain[3][ele.idx[j]] + lambda *
			// (strain[0][ele.idx[j]] + strain[1][ele.idx[j]]);
			// } else {
			// tx = (lambda + 2 * mu) * strain[0][ele.idx[j]] + lambda *
			// strain[1][ele.idx[j]];
			// ty = lambda * strain[0][ele.idx[j]] + (lambda + 2 * mu) *
			// strain[1][ele.idx[j]];
			// stress[0][ele.idx[j]] += tx;
			// stress[1][ele.idx[j]] += ty;
			// if (BVPType.PLANESTRAIN == tp) {
			// stress[3][ele.idx[j]] += Nu * (tx + ty);
			// } else if (BVPType.PLANESTRESS == tp) {
			// strain[3][ele.idx[j]] += -Nu * (tx + ty) / E;
			// }
			// }
			// stress[2][ele.idx[j]] += mu * strain[2][ele.idx[j]];
			// angleCount[ele.idx[j]] += angle[j];
			// }
			// }
			// for(int i = 0; i < n; i++) {
			// strain[0][i] /= angleCount[i];
			// strain[1][i] /= angleCount[i];
			// strain[2][i] /= angleCount[i];
			// strain[3][i] /= angleCount[i];
			// stress[0][i] /= angleCount[i];
			// stress[1][i] /= angleCount[i];
			// stress[2][i] /= angleCount[i];
			// stress[3][i] /= angleCount[i];
			// }

			for (int j = 0; j < 3; j++) {
				strain[0][ele[j]] += dxU * angle[j];
				strain[1][ele[j]] += dyV * angle[j];
				strain[2][ele[j]] += (dxV + dyU) * angle[j];
				if (BVPType.AXISYMMETRIC == tp) {
					strain[3][ele[j]] += eleU[j] / x[j] * angle[j];
				}
				angleCount[ele[j]] += angle[j];
			}
		}
		for (int i = 0; i < n; i++) {
			strain[0][i] /= angleCount[i];
			strain[1][i] /= angleCount[i];
			strain[2][i] /= angleCount[i];
			strain[3][i] /= angleCount[i];
		}

		for (int i = 0; i < n; i++) {
			if (BVPType.AXISYMMETRIC == tp) {
				stress[0][i] = (lambda + 2 * mu) * strain[0][i] + lambda * (strain[1][i] + strain[3][i]);
				stress[1][i] = lambda * (strain[0][i] + strain[3][i]) + (lambda + 2 * mu) * strain[1][i];
				stress[3][i] = (lambda + 2 * mu) * strain[3][i] + lambda * (strain[0][i] + strain[1][i]);
			} else {
				stress[0][i] = tx = (lambda + 2 * mu) * strain[0][i] + lambda * strain[1][i];
				stress[1][i] = ty = lambda * strain[0][i] + (lambda + 2 * mu) * strain[1][i];
				if (BVPType.PLANESTRAIN == tp) {
					stress[3][i] = Nu * (tx + ty);
				} else if (BVPType.PLANESTRESS == tp) {
					strain[3][i] = -Nu * (tx + ty) / E;
				}
			}
			stress[2][i] = mu * strain[2][i];
			meanStress = (stress[0][i] + stress[1][i] + stress[3][i]) / 3;
			t1 = stress[0][i] - meanStress;
			t2 = stress[1][i] - meanStress;
			t3 = stress[3][i] - meanStress;
			misesStress[i] = 0.5 * (t1 * t1 + t2 * t2 + t3 * t3) + stress[2][i] * stress[2][i];
			misesStress[i] = Math.sqrt(3.0 * misesStress[i]);
		}
	}

}
